clc
clear

ncfile = '../run/hope_output.n90.nc';

var_name  = 'z';
LevelStep = 40;
cmin = 4800;
cmax = 6000;

% var_name  = 'rho';
% LevelStep = 0.0025;
% cmin = 0.5;
% cmax = 1.5;

% var_name  = 't';
% LevelStep = 1;
% cmin = 232;
% cmax = 272;

% var_name  = 'pt';
% LevelStep = 2;
% cmin = 280;
% cmax = 332;

% var_name  = 'mr';
% LevelStep = 0.001;
% cmin = 0;
% cmax = 0.02;

% var_name  = 'us';
% LevelStep = 4;
% cmin = -20;
% cmax = 20;

% var_name  = 'vs';
% LevelStep = 4;
% cmin = -20;
% cmax = 20;

% var_name  = 'w';
% LevelStep = 0.02;
% cmin = -1;
% cmax = 1;

% var_name  = 'p';
% LevelStep = 50;
% cmin = 4800;
% cmax = 6000;

plt_res   = ncreadatt(ncfile,'/','dx') / 2;
plev      = 50000; % select pressure level to plot, in Pa
ilev      = 1; % For plot pressure only
FontSize  = 6;
it        = 2;

lon = ncread(ncfile,'lon');
lat = ncread(ncfile,'lat');
p   = ncread(ncfile,'p',[1,1,1,1,it],[Inf,Inf,Inf,Inf,1]);
if(~strcmp(var_name,'z'))
    var = ncread(ncfile,var_name,[1,1,1,1,it],[Inf,Inf,Inf,Inf,1]);
else
    var = ncread(ncfile,var_name,[1,1,1,1,1],[Inf,Inf,Inf,Inf,Inf]);
end
    
lon_min = 0;
lon_max = 360;
lat_min = -90;
lat_max = 90;

nx = size(lon,1);
ny = size(lon,1);
np = size(lon,3);

nz = size(var,4);
nt = size(var,1);

var_p = zeros(nx,ny,np); % variable on pressure level
for ip = 1:np
    for j = 1:ny
        for i = 1:nx
            if(strcmp(var_name,'z'))
            var_p(i,j,ip) = interp1(squeeze(log(p(i,j,:,ip))),...
                                    squeeze(var(1,i,j,:,ip)),log(plev));
%             var_p(i,j,ip) = spline(squeeze(log(p(i,j,:,ip,it))),...
%                                     squeeze(var(1,i,j,:,ip)),log(plev));
%             var_p(i,j,ip) = interp1(squeeze(p(i,j,:,ip,it)),...
%                                     squeeze(var(1,i,j,:,ip)),plev);
            elseif(~strcmp(var_name,'p'))
              var_p(i,j,ip) = interp1(squeeze(log(p(i,j,:,ip))),...
                                      squeeze(var(i,j,:,ip)),log(plev));
%               var_p(i,j,ip) = spline(squeeze(log(p(i,j,:,ip,it))),...
%                                       squeeze(var(i,j,:,ip,it)),log(plev));
%               var_p(i,j,ip) = interp1(squeeze(p(i,j,:,ip,it)),...
%                                       squeeze(var(i,j,:,ip,it)),plev);
            end
        end
    end
end
if strcmp(var_name,'p')
    var_p = squeeze(var(:,:,ilev,:));
end

lon1d = reshape(lon,[],1);
lat1d = reshape(lat,[],1);

lon_plt = lon_min:plt_res:lon_max;
lat_plt = lat_min:plt_res:lat_max;
[lon_plt,lat_plt] = meshgrid(lon_plt,lat_plt);

var1d = reshape(var_p,[],1);

% extend for interpolation
n1d = size(lon1d,1);
lon1d_ext(0*n1d+1:1*n1d) = lon1d - 360;
lon1d_ext(1*n1d+1:2*n1d) = lon1d;
lon1d_ext(2*n1d+1:3*n1d) = lon1d + 360;
lat1d_ext(0*n1d+1:1*n1d) = lat1d;
lat1d_ext(1*n1d+1:2*n1d) = lat1d;
lat1d_ext(2*n1d+1:3*n1d) = lat1d;
var1d_ext(0*n1d+1:1*n1d) = var1d;
var1d_ext(1*n1d+1:2*n1d) = var1d;
var1d_ext(2*n1d+1:3*n1d) = var1d;

var_plt = griddata(lon1d_ext,lat1d_ext,var1d_ext,lon_plt,lat_plt,'linear');
% var_plt = griddata(lon1d_ext,lat1d_ext,var1d_ext,lon_plt,lat_plt,'cubic');

figure
% % pcolor plot
% plt = pcolor(lon_plt,lat_plt,squeeze(var_plt));
% set(plt,'EdgeColor','None')
% % shading interp
% % set(gca,'Clim',[cmin,cmax])
% colormap(jet)
% colorbar

% contour plot
[plt,h]=contour(lon_plt,lat_plt,squeeze(var_plt),'LevelStep',LevelStep,'ShowText','on');
% set(gca,'Clim',[cmin,cmax])
clabel(plt,h,'FontSize',FontSize)
colormap(jet)

title([var_name,' at ',num2str(plev),'Pa'])
varname_file = strrep(var_name,'\','');
% print(gcf,'-r600','-dpng',['HOPE_',varname_file,'_at_',num2str(plev),'Pa','.png']);